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Abstract. We develop a formalism for the calculation of the macroscopic 
dielectric response of composite systems made of particles of one material 
embedded periodically within a matrix of another material, each of which 
is characterized by a well defined dielectric function. The nature of these 
dielectric functions is arbitrary, and could correspond to dielectric or conducting, 
transparent or opaque, absorptive and dispersive materials. The geometry of 
the particles and the Bravais lattice of the composite are also arbitrary. Our 
formalism goes beyond the longwavelenght approximation as it fully incorporates 
retardation effects. We test our formalism through the study the propagation 
of electromagnetic waves in 2D photonic crystals made of periodic arrays of 
cylindrical holes in a dispersionless dielectric host. Our macroscopic theory 
yields a spatially dispersive macroscopic response which allows the calculation of 
the full photonic band structure of the system, as well as the characterization 
of its normal modes, upon substitution into the macroscopic field equations. 
We account approximately for the spatial dispersion through a local magnetic 
permeability and we analyze the resulting dispersion relation, obtaining a region 
of left- handedness. 



1. Introduction 

The propagation of light within homogeneous materials is completely characterized by 
their electromagnetic response. In the optical regime, magnetic effects are typically 
negligible so that knowledge of the dielectric response is sufficient for the study of 
wave propagation L. However, when the system has spatial inhomogeneities, the 
behaviour of light is not trivial and has captured the attention of many researchers 
[21E1I11E1IS1EIIH1E1 [IHl HI]- The possibility of controlling the propagation of photons 
by producing artificially structured materials has been widely demonstrated in recent 
years. Novel effects have been obtained, such as negative refraction, inverse Doppler 
effect, optical invisibility cloaking, optical magnetism etc. [I2l [13l [HI [15] . 

Effective medium theories [El[ni[H[H[2Ql[2Tl[22l[23l[24] have been 
proposed to describe inhomogeneous systems, such as diluted colloidal suspensions 
in the quasistatic or long wavelength limit, in terms of a homogeneous macroscopic 
response. Further developments in homogenization of composites have been proposed 
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[231|2Bl|271|2Sl|5Sl|3Dl|3Tll31l33]and their limits of validity have have been discussed 
[3 GH |35]. 



In this paper we obtain the macroscopic optical response using a computationally 
efficient reformulation of a procedure [SJ |5] original developed to account for the 
local field effect of systems with spatial fluctuations. Our main result is that for 
any system, the macroscopic response to an external excitation is simply the average 
projection of the corresponding microscopic response. We apply this general result 
to Maxwell equations in order to obtain an explicit expression for the macroscopic 
dielectric response of binary composites. We illustrate its use with numerical results 
for a 2D periodic photonic crystal. As the wavelength of the electromagnetic field 
becomes comparable to the characteristic microscopic spatial lengthscale of the system, 
the macroscopic response becomes highly non-local, with a strong dependence on the 
wavevector besides its usual dependence on frequency. When the full spatial dispersion 
is taken into account, the photonic band structure of the system can be obtained and 
analyzed from its macroscopic response. 

The paper is organized as follows: In Sec. [2] we formulate a general theory for 
the macroscopic dielectric response of composite systems. In Sec. |3j we adapt our 
formulation to periodic systems made of two alternating materials, each characterized 
by well defined dielectric functions. In Sec. |4] we develop an efficient computational 
scheme which allows the numerical calculation of the macroscopic response without 
recourse to explicit operations between large matrices. In Sec. [5jwe obtain the non- 
local macroscopic dielectric tensor and use it to calculate macroscopically the full band 
structure of a 2D photonic crystal made up of non-dispersive dielectric components; 
we compare it to exact results as well as to results within a local approximation that 
partially accounts for spatial dispersion through an effective magnetic permeability. 
Finally, in Sec. [6] we present our conclusions. 

2. General theory 

Consider a non-magnetic inhomogeneous system characterized by its dielectric 
response e, defined through the constitutive equation D = eE, where D and E are the 
microscopic displacement and electric fields respectively. We designate these fields 
as microscopic as they have spatial fluctuations due to the inhomogeneities of the 
material. They are not microscopic in the atomic scale, but rather, on the scale of 
the spatial inhomogeneities of the system. We use a caret (") over a symbol to denote 
its operator nature and we leave implicit the dependence of the dielectric response on 
position as well as on frequency. Our purpose is to obtain the macroscopic dielectric 
response of the system, relating the macroscopic displacement and electric fields, from 
which the fluctuations have been removed. 

We start from Maxwell equations for monochromatic microscopic fields. We follow 
the usual procedure [33] to decouple the magnetic field from the electric field to obtain 
the second order wave equation 



WE 



4ttc 
iq 
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is the wave operator, q — uj/c is the wavenumber of light in vacuum, uj is the frequency, 
c is the speed of light in vacuum and J is the external electric current density. Here we 
introduced the transverse projector V T such that the transverse projection of a field 
F is F T = V T F. For completeness, we also introduce the longitudinal projector V L 
such that V L + V T = 1, with 1 the identity operator. Notice that Eq. (fTJ) contains 
both a longitudinal and a transverse part, so it describes not only transverse waves, 
but also allows for the possible excitation of longitudinal waves such as plasmons [37 . 
We solve Eq. ((T|) formally for E to obtain 

47TC - -, 

E = W _1 J. (3) 

iq 

As we are interested on the macroscopic response of the system, we introduce the 
average V a and fluctuation Vf projectors, such that an arbitrary field F can be 
written as F = F Q + F/ with F a = V a F its average and F/ = VfF its fluctuating 
part, and we identify the average field with the macroscopic field F a = Fjf. Later, 
we will provide appropriate explicit definitions for V a and Vf] here we remark that 
they must be idempotent V 2 — T^a, Vf — Vf, i.e., the average of the average is the 
average, and the fluctuations of the fluctuations are the fluctuations 5 . Furthermore, 

PaVf = VfVa = and V a +Vf = l. 

As the macroscopic response would be useless unless the external excitations are 
devoid of microscopic spatial fluctuations, we assume that the external current has no 
fluctuations, so J/ = and J = J M = V a 3 M , where we used the idempotency of V a . 
Thus, acting on both sides of Eq. ([3]) with V a we obtain 

EM = ^- JM (4) 

Here we define O a p = V a OVfs, with a, f3 = a, / for any operator O. As Eq. (UJ) 
relates the macroscopic external electric current to the macroscopic electric field we 
may identify the macroscopic inverse wave operator Wj^ 1 , 

iq 

where 

Wm^W- 1 . (6) 

We summarize this results stating that the macroscopic inverse wave operator is simply 
the average of the microscopic inverse wave operator. This is a particular case of a 
more general result: the macroscopic response to an external excitation is simply the 
average of the corresponding microscopic response. 

From the macroscopic Maxwell equations we may further relate the macroscopic 
wave operator in Eqs. ([5]) and (J6j) to the macroscopic dielectric response of the system 
e M through 

f 

Q 

in analogy to Eq. ([2]). Thus, we have to invert the wave operator, average it, and 
invert it again to finally identify the macroscopic dielectric response. Notice that we 
have made no approximation whatsoever. We remark that as e relates two fields, 
E and D, that have spatial fluctuations, we may not simply average e to obtain its 
macroscopic counterpart, i.e. e M ^ e aa . The difference constitutes the local field effect 
0. 



W M = e M + —V 2 VaP T , (7) 
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Our result may easily be generalizable to other situations and other response 
functions. The procedure consists on first identifying the response (in our case VV _1 ) 
to the external perturbation (i.e. J), and then averaging it to yield the corresponding 
macroscopic response (i.e. W^ 1 = VV^ 1 ), which may further be related to the desired 
macroscopic response operator (i.e. e M ). 



3. Periodic binary systems 

In this section we use Eq. ((6]) to obtain the optical properties of an artificial binary 
crystal made of two materials A and B with dielectric functions £a and es . We assume 
that both media are local and isotropic so that ca and es are simply complex functions 
of the frequency. For convenience, we will further assume that is real, though this 
assumption is easily relaxed |38j . 

We introduce the characteristic function B(r) of the inclusions, such that B(v) = 1 
whenever r is on the region B occupied by the inclusions, and B(r) = otherwise. 
Thus, we may write the microscopic dielectric response as 

e(r) = H( u _fl( r )), (8) 
u 

where we defined the spectral variable u = 1/(1 — cb/^a) [17] . The microscopic wave 
operator of Eq. @ may be written as 



W 



(ug- 1 - B) , (9) 

where the characteristic operator B corresponds to multiplication by B(r) in real 
space, and we defined 

■ <io) 

which, as shown below, plays the role of a metric. 
Using Eq. © we write Eq. ([6]) as 

Wm 1 = —9aa(u-BgY\ (11) 

V / aa 

where we have taken advantage of the fact that g is unrelated to the texture of the 
crystal, so that it does not couple average to fluctuating fields. 

Using Bloch's theorem [35] > we consider an electric field of the form [TO] 

E k (r) -^E G e J ( k+G > r , (12) 

G 

where k is a given wavevector, {G} is the reciprocal lattice of our crystal and the 
coefficients Eq represent the field in reciprocal space. In this representation, all 
operators may be written as matrices with vector index pairs G,G', besides other 
possible indices, such as Cartesian ones. Thus, we represent the longitudinal projector 
as the matrix 

L (k + G) (k + GQ 

|k + G| |k + G'| (13) 

with <5gc the Kronecker's delta, so the transverse projector becomes V^ G , = 
l^GG' ~ Pqq; with 1 the Cartesian identity matrix. The Laplacian operator in 
reciprocal space is represented by 

V 2 -> -(k + G) 2 ,W, (14) 
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and we can define the average projector as a cutoff in the reciprocal space [B] 

(Po)gg- = <5go<5g'o, (15) 

so that average fields simply keep the contributions with vector k while all other 
wavevectors are filtered out. 

As shown by Eq. dill) , we only require 

{u-Bg)-i = {u-Bg)zl (16) 

to obtain the macroscopic inverse wave operator, where the subindices denote the 
projection onto the subspace with G = 0. 

4. Recursive method 



The calculation of Eq. ([16]) is analogous to that of a projected Green's function |40j |4T] 
Gaa(e) = (a\G(e)\a), (17) 

onto a given state \a), where 

g{e) = {e-H)- 1 (18) 

is the Green's operator corresponding to some Hermitian Hamiltonian % and e is a 
complex energy. In Ref. |26j a similar result was obtained, where the Hamiltonian H 
was identified with the longitudinal projection of the characteristic function B LL , the 
energy s with the spectral variable it, and |a) with a slowly varying longitudinal wave, 
and it was shown that the corresponding projected Green's function was proportional 
to the inverse of the longitudinal macroscopic dielectric function. Haydock's method 
may be applied to obtain the projected Green's function (|17j) in a very efficient way 
[42] [43] . and it has been adapted previously [26] to the calculation of the optical 
response of nanostructured systems in the long-wavelength local limit. In this section 
we generalize the approach of Ref. |26] to arbitrary frequencies and wavevectors. 

Eqs. (fT7]) and (fT8]l arc similar to Eq. (fTo) , although the operator product 
H = Bg that plays the role of Hamiltonian does not correspond to a symmetric 
matrix. Nevertheless, it would correspond to a Hermitean operator if we use g as 
a metric, that is, if we use {tp\g\(f>) instead of (iplfi) as the scalar product of two 
arbitrary states \tp) and \<j>). Then, it is easy to verify the Hermiticity condition 
KV-'I'K'^I^))]* = i^Qv^W))- Notice, however, that g is not positive definite. 

According to Eq. ([TT]) . we need the average (u — Bg)^ — > (a\(u — i3g) _1 |a), where 
we projected onto an average state \a) — &o|0) consisting of a plane wave with a given 
wavevector k, frequency uj and polarization e. Here, |0) is normalized according to 
the metric g, i.e., (0|.g|0) = go = ±1, and the coefficient &o is a real positive number 
chosen to guarantee that \a) is normalized in the usual sense (a\a) = 1. Now, we 
define | — 1) = and we obtain new states by means of the recursion relation 

\n+ 1) = Bg\n) = b n+1 \n + 1) + a n \n) + g n ~ig n b n \n - 1), (19) 

where all the states \n) are orthonormalizcd according to the metric <?, that is 

(n\g\m) = g n S„ m , (20) 

with g n = ±1 and 5 nm the Kronecker's delta function. The requirement of 
orthonormality yields the generalized Haydock coefficients a n , b n +i, and g n +i given 
the previous coefficients b n , g n and g n -i- Thus, a n are obtained from 

(n\g\n + 1) = a n g n , (21) 
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and b n+ i from 

' ' ' ~~ " ' J n+1 



l\g\n + 1) = g n +ib 2 n+1 + g n a 2 n + g„-ib 2 n , (22) 



where we choose the sign g n +i — ±1 so that b 2 l+1 is positive and we may choose 6„+i 
as a real positive number. In the basis {!«)}, Bg is represented by a tridiagonal matrix 
with a n along the main diagonal, b n along the subdiagonal and g n -ignb n along the 
supradiagonal, so that 



Bg^ 



I u- a -hgxgo 

-&i u-ai -b 2 g 2 gi 

-62 u - a 2 -b 3 g 3 g2 

V ; ; ; 



\ 



(23) 



Notice that we may represent the states |n) through the corresponding Cartesian 
field components for each reciprocal vector G or for each position r within a unit 
cell. Thus, the action of g is a trivial multiplication in reciprocal space, while that 
of B is a trivial multiplication in real space, and we may repeatedly compute Bg\n) 
and Haydock's coefficients without performing any actual matrix multiplication, by 
alternatingly fast-Fourier transforming our representation between real and reciprocal 
space. 

According to Eq. (fTTj) . we do not require the full inverse of the matrix (1231) . but 
only the element in the first row and first column. Following Ref. |26j , we obtain that 
clement as a continued fraction, which substituted into Eq. (JTTJ) yields 

a A u — »o ^2 

u—ai — -5- 

„_a 2 i 

The right hand side of Eq. (j2"4")l denotes the macroscopic response projected onto 
a state with the given wavevector k, frequency u and polarization e. Due to the 
translational invariance of the homogenized system, the inverse wave operator for a 
given k corresponds simply to a tensor with Cartesian components {WZr)ir Thus, 
so far we have calculated its inner products with e, as shown by the second term 
of Eq. (|2"4")l . We may repeat the calculation above for different orientations of the 
(possibly complex) unit polarization vector e and solve the resulting equations for all 
the individual Cartesian components (W^ 1 )^-. Finally, we perform a matrix inversion 
of this macroscopic tensor and we use Eq. [7] and [5] to compute the macroscopic 
dielectric tensor 

<#(w, k) = l(fc 2 ^- - hk 3 ) + k). (25) 

Eqs. (|24|) and (f25"|) constitute the main result of our formalism. Notice that in general, 
this dielectric tensor would depend on both uj and k, that is, it is a temporal and 
spatially dispersive response function. 

5. Results 

To test our formalism, we first calculate the local limit of the longitudinal macroscopic 
dielectric function e^ 1 (uj) = e*^(w, k — > Ox) of a 2D system made up of a square array 
of empty cylindrical holes (e# = 1) placed with lattice parameter a within a dielectric 
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Figure 1. (Color online) Longitudinal macroscopic dielectric function (ui) 
of a 2D system (inset) made up of a square array with lattice parameter a of 
empty cylindrical holes (eg = 1) of circular cross section with radius p = 0.45a 
within a dielectric medium with permittivity €a = 12. We show the results of 
using Maxwell-Garnett formula (MG), of a nonretarded calculation (NR), of a 
numerically exact matrix calculation (diamonds) and of our theory (solid lines) 
for two different lattice parameters a = 40nm and 120nm. 



medium with permittivity = 12. The holes are of circular cross section with 
radius p = 0.45a. We assume that the axes of the cylinders are parallel to the z 
axis. This system has been frequently been used as a testground for calculations of 
photonic crystal properties [44.. In this calculation we took the limit k > along 
the x direction, and we introduced the Cartesian unit vectors x, y and z. We remark 
that it is important to assign a direction to k even in this limit, as the transverse T 
and longitudinal L projection of W M differ, although, for this system, the resulting 
e M (ui) = e M (w,k — > 0) is isotropic within the xy plane. 

In Fig. [T]we show e^ 1 (uj) calculated for two different lattice parameters a = 40nm 
and a = 120nm. Our calculations were performed using the Perl Data Language 
[3S] U . The figure shows that our formalism yields the same results as a full 
straightforward matricial calculation [3D] that solves Maxwell equations in a plane 
wave basis. Nevertheless, our calculation is much more economical in memory usage, 
as we don't have to store the full matrices that represent the dynamics of the system, 
and is at least four orders of magnitude faster as we do not perform any matrix 
manipulation with our scheme. To illustrate the effects of retardation, in the figure we 
also show the result of using the non-retarded calculation of Ref. [26] (NR), together 
with those of the well-know 2D Maxwell-Garnet (MG) effective medium formula [19 . 
As the dielectric functions of the components are independent of the frequency, since 
the constituents were taken to be non-dispersive, in this case the frequency dependence 
of the result arises solely from to the finite ratio of the free space wavelength to the 
lengthscale of the system, namely, the lattice parameter a. Thus, in the u) — > limit 
our results approach those of the NR calculation. The difference between the NR 
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Figure 2. (Color online) Out-of-plane macroscopic dielectric function ejf (a;, k/\) 
of the same system as in Fig. [T]for a given wavevector along the A line of the 
first Brillouin zone (solid lines). The line (k&/q) 2 (long dashes) intersects (circles) 
e^f (id, k/\) at the frequencies uj a (]<./\) of the TE normal modes, with a =1, II, IV, 
V. The intersections (diamonds) of e^f(id,k^) (short dashes) with (k'^/q) 2 (long 
dashes) corresponding to a shifted wavevector k^ = k/\ + (27r/a)y yields the 
modes a =111, IV (see text). 



and MG calculations is due to the high filling fraction / w 0.64 of the cylinders, 
which interact non-dipolarly with their nearest neighbors which they almost touch, 
while MG contains only dipolar interactions. Naturally, retardation effects become 
stronger as a increases. Furthermore, e M (w,k) depends only on the products qa 
and ka so that the two curves shown in Fig. [T] could be collapsed into one curve 
if appropriately plotted. We use this in Fig. [2j which shows the transverse response 
elf 1 (w, kA) = (u>, kA) corresponding to out-of-plane or TE polarization as a function 
of the normalized frequency qajlix. Here, we have chosen a finite in-plane wavevector 
kA = ( 7r /2a,0, 0), corresponding to the midpoint of the A line between V and X in 
the first Brillouin zone (BZ) [35] ■ This figure covers a larger frequency range than 
Fig. [TJ and therefore it displays a series of poles related to resonant multiple coherent 
reflections at the interfaces between the A and B regions. 

The dispersion relation of transverse waves with TE polarization propagating 
along the xy plane can be simply written as |46j 

4V,k) = ^. (26) 
T 

Therefore, in Fig. [5]we also plotted the curve (k&/q) 2 . Those points where it intersects 
the curves for correspond to normal TE modes of the system with frequencies 
uj a (k = kA), a —I, II, ...Here, we use roman numerals to number the modes in 
increasing frequency order. 

In Fig. [3J we show the photonic bands (circles) of the TE modes of our 
system obtained by solving the dispersion relation ([2"B"f as k travels along the line 
r — A — X — Z — M — £ — T within the first BZ. In order to test our theory, we also 
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Figure 3. (Color online) Photonic bands of the TE modes of the same system 
as in Fig. [l] obtained by solving Eq. I|26ll (circles). For comparison, we also show 
the photonic bands obtained from the solution of the corresponding eigenvalue 
problem, Eq. H27H , as described in the text (solid lines). We also show the modes 
obtained by shifting the wavevector k — > k' = k + (2n/a)y (diamonds). 



plot the TE modes obtained by solving the eigenvalue equation [47] 




through an implicitly restarted Arnoldi iteration [48| , applying 1 / <\J e(r) in real space 
and V — >• i(k + G) in reciprocal space [38]. We further verified our results by 
recalculating them with the freely available package MPB [50] and comparing it to 
previous results [44] for the same system. The agreement between these calculations 
shows that it is feasible to use the macroscopic dielectric response of the system to 
obtain its photonic band structure. However, to succeed, we have to account fully for 
the spatial dispersion of e^f(uj, k). 

We remark that although we obtained the correct photonic bands even for large 
wavevectors, going all the way around the first BZ, there are some regions where we 
failed to obtain all of the normal modes. For example, our procedure did not produce 
the third band (labeled III in the figure) along the A line, and neither the fourth 
band (IV) along the E line, between M and T. Furthermore, it did not yield the 
second band (II) which is almost degenerate with the third band (III) along E. It is 
easily shown that the microscopic field corresponding to the missing band along A is 
antisymmetric with respect to the reflection G y <H> — G y about the A line. Similarly, 
the missing bands along the region E are antisymmetric with respect to the reflection 
G x «-» G y about the E line. Thus, the reciprocal vector does not contribute to the 
electric field for those bands [51], i.e., the missing modes have no macroscopic field 
components and thus, apparently they may not be obtained from the roots of the 
macroscopic dispersion relation (|26|) . 

Nevertheless, within our formalism, the wavevector k is actually not restricted to 
lie within the first BZ. Thus, we may calculate the macroscopic dielectric function for 
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wavevectors beyond the first BZ, and obtain the modes in the extended scheme. In Fig. 
Owe show part of the macroscopic dielectric function e^(uj,k' A ) with a wavevector 
that has been shifted out of the first BZ from along the y direction by the reciprocal 
vector (27r/o)y, i.e., k' A = kA + (2n/a)y. We remark that e M (w,k) is not a periodic 
function of k, as it corresponds to a specific plane wave, not to a microscopic Bloch's 
wave. The intersection of the curve (k' A /q) 2 with e^(uj,k' A ) yields the corresponding 
macroscopic modes. For example, the diamonds in Fig. [2] illustrate two of the normal 
modes that could be obtained using this shifting procedure. One of these modes is 
identical to the mode labelled IV that was obtained previously using the unshifted 
response '(u) ,\ca) • Nevertheless, another mode, labelled III appears just below 
qa /2tt = 0.4 and it corresponds to one of the previously missing modes. We may 
shift it back into the first BZ in order to display it in the reduced zone scheme. 
Proceeding in this fashion, we have obtained all of the previously missing bands, 
shown by diamonds in Fig. [3J Thus, we have shown that we can calculate the full 
photonic band structure for the TE modes from the macroscopic non-local dielectric 
function of the composite system. 

A similar procedure to that discussed above can be employed also for the TM 
modes of the system, with the electric field within the xy plane. In Fig. 0]we show the 
in-plane macroscopic dielectric functions c^,(u),\s.a) and Cyy(uj,k A ) of the system for 
the same wavevector kA as in Fig. [5] This wavevector and the system are symmetric 
under y <-> —y reflections. Therefore, e^f y = and we may identify the transverse 
and longitudinal responses as e^(uj,k A ) = e^(w, kA) and ej^(uj,k A ) = €^.(u),\s.a) 
respectively. The normal modes are then obtained from the transverse dispersion 
relation (f26| and the longitudinal dispersion relation [46] 

e£V,k)=0. (28) 

These are identified with circles and diamonds in the figure. 

In Fig. [5] we do a similar calculation but along the £ line, for the wavevector 
ks = (7r/2a, n/2a). Along this line, the wavevector and the system are symmetric 
under the reflection x y, and thus e^(w,ks) = e^(w,ks). We may now identify 
the transverse and longitudinal responses as e^(ui, kg) = e*^(w, ks) — e^(w, kg) and 
e^(oj,ks) = e^(w, kg) +e^,(u;,ks) respectively. The upper panel of Fig. [5] shows 
£y f (cj,ks), (ks/q) 2 and their intersections (circles) which yield the transverse modes, 
while the lower panel displays e^ 1 (oj, ks) and its zero (diamond) which correspond to 
a longitudinal mode. 

Along the Z line, from X to M , there is no crystal symmetry operation that 
leaves invariant the wavevector, and thus transverse and longitudinal fields mix among 
themselves, i.e., there are no longitudinal nor transverse modes. Nevertheless, we 
may obtain the frequencies of the mixed modes through the singularities of the wave 
operator matrix 

det [W M (uj,k)] = 0, (29) 

as illustrated in Fig. |6]for a wavevector kz = (tt / 'a, 7r /2a) along the Z line. 

From calculations such as those illustrated by Figs. [U [5] and El or rnore generally, 
from the zeroes of det [W M (w, k)] for arbitrary wavevectors k, we can obtain the full 
TM photonic band structure, as shown by Fig. [7] In order to test our theory, we also 
plot the TM modes obtained by solving the eigenvalue equation [47] 

V'^VB,(r) = -^B,(r) (30) 



Macroscopic optical response and photonic bands 

14 
12 



11 



< 



10 

8 
6 
4 
2 

-2 
-4 



1 1 
- \ (k A /q) 2 


ii i 




\ I ^ 








II / 


IV/ 


i i ifi A 



14 
12 
10 
8 
6 
4 
2 

-2 
-4 



1 1 1 1 


1 1 

III// 




1 1 1 1 / 







0.1 



0.2 



0.3 

qa/2n 



0.4 0.5 0.6 



Figure 4. (Color online) In-plane macroscopic dielectric functions eif(u), Ica) 
(upper panel) and eJ^(uj,kA) (lower panel) of the same system as in Fig. [T]for 
the same wavevector as in Fig. [2] The line (k&/q) 2 (dashes) intersects (circles) 
e^f(o;,kA) at the frequencies aj a (kA) of the TM transverse normal modes, with 
a =1, II, IV. The zeroes (diamonds) of e^(oj, kA) yield the longitudinal mode III 
(see text). 



using similar methods [48, ^38] as for the TE case and comparing them successfully to 
previous results [44] for the same system. The agreement between these calculations 
shows that it is also feasible to use the macroscopic dielectric response of the system 
to obtain its photonic band structure in the TM case. 

An approximate band structure could be obtained for transverse waves in the 
region of small wavevectors k — > by expanding the LHS of Eq. (f26|) up to second 
order in k and solving for k 2 . The result is a local dispersion relation 

k 2 = q 2 e M (oj)fi(Lo) 1 (31) 

where the non- locality of e^ f (u;,k) is partially accounted for through an effective 
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Figure 5. (Color online) In-plane macroscopic dielectric functions eir (a;, kj;) 
(upper panel) and e* f (cj,kx;) (lower panel) of the same system as in Fig. [T] 
for a given wavevector kj; along the £ line. The line (k-^/q) 2 (dashes) intersects 
(circles) e^f (ui, ks) at the frequencies u) a (k^) of the TM transverse normal modes, 
with a=I, II, IV. The zeroes (diamonds) of (u), ks) yield the longitudinal mode 
III (see text). 



magnetic permeability 



1 



1 



(32) 

2 dk 2€ T (W,KXJ fe=Q 

In Fig. [8] we show the resulting approximate bands for the TE case. We can see that 
the acoustic band within the local approximation is indistinguishable from the exact 
result for small wavevectors. This local acoustic band extends up to the frequency 
qa/2ir w 0.2 for which [i has a pole. The local approximation also reproduces the upper 
band of the exact calculation (labelled V in Fig. \S§ , but only in the immediate vicinity 



of T, close to the zero qa/2ir ss 0.43 of e 



M, 



For larger wavevectors it acquires a 
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Figure 6. Determinant of the macroscopic wave operator matrix W M (cj,kz) 
for a given wavevector along the Z line. The normal modes of the system 
(iia(kz), a =1, II, III, IV, correspond to its zeroes (squares). 




Figure 7. (Color online) Photonic bands of the TM modes of the same system as 
in Fig. [T] The modes between Y and X and between M and Y may be obtained 
by solving Eq. I|26j l for the transverse (T) modes (circles) and Eq. I|28|l for the 
longitudinal (L) modes (diamonds), or Eq. (1 29 I t for both of them. The modes 
from X — > M are obtained from Eq. H29H and have a mixed polarization (squares). 
For comparison, we also show the photonic bands obtained from the solution of 
the corresponding eigenvalue problem, Eq. H30I I. as described in the text (solid 
lines) . 
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Figure 8. Local dispersion relation for the same system as in Fig. [T]for the TE 
case (solid lines, left panel), obtained from Eq. (13 1 D - Local dielectric function 
e M (uj) (middle panel). Local magnetic permeability (solid line, right panel) 

and 10 4 /i(o;) (dash-dotted line, right panel). For comparison purposes, we also 
show the exact dispersion relation (dashed line, left panel). 



larger positive dispersion. As discussed above, the band labelled III in Fig. [3] doesn't 
couple to macroscopic fields for k along the A line, while the band labelled IV doesn't 
couple along the E line. Thus, neither band couples to macroscopic fields at T and as 
a consequence, these bands are not reproduced at all by the local approximation (|3lj). 

In the region qa/2n « 0.36 — 0.40 the local approximation predicts a curious 
band that displays backbending. Its lower part begins at the pole of e M (w), which 
corresponds to a double zero of //(w) [38]. This part of the local dispersion relation 
is spurious, as the second order Taylor expansion of the dielectric function that 
yields Eqs. (|32|) and (|3"Tj) starting from Eq. (|26|) would be meaningless at a pole. 
Consequently, this part of the local band doesn't correspond to any band within the 
exact calculation. On the other hand, the upper part of the backbending band starts 
at a simple zero of p, which arises from a different kind of singular behavior, consisting 
of a pole in e|/(w,k) at qa/2ir « 0.40 which is suppressed at k — as its weight is 
approximately proportional to fc 2 . This part of the local dispersion relation agrees with 
band II of Fig. [3]at T and displays a negative dispersion as does band II. Furthermore, 
for this band, both e M and /i are negative. Thus, our photonic crystal, made up of 
holes within a dispersionless dielectric, behaves for some frequencies as a left-handed 
metamaterial [52l [53] . Nevertheless, caution should be exercised when using Eqs. (|32|) 
and (|3"Tj) close to singularities of e^(uj, k). 

6. Conclusions 



We have shown that the macroscopic inverse wave operator of composite systems 
with arbitrary spatial fluctuations is simply given by the average projection of 
the corresponding microscopic operator. These operators may be related to the 
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macroscopic and macroscopic dielectric functions, thus yielding a general procedure to 
incorporate the effects of spatial fluctuations in the calculation of all the components of 
the macroscopic dielectric tensor of the system. We have extended Haydock's recursive 
scheme in order to calculate very efficiently the response of periodic two-component 
systems in terms of a continued fraction whose coefficients may be obtained without 
recourse to operations with the large matrices that characterize the microscopic 
response and the field equations. We illustrated and tested our results through 
the calculation of the wavevector and frequency dependent macroscopic dielectric 
tensor of a 2D dielectric system with non-dispersive non-dissipative components. We 
identified the transverse and longitudinal components of the response for wavevectors 
along symmetry lines and we obtained all the components of the tensorial response 
for the case of reduced symmetry. The macroscopic response has a series of poles, 
related to resonant multiple coherent reflections and by accounting for its spatial 
dispersion we obtained the full photonic band structure of the system from the 
dispersion relations corresponding to the homogenized material. Besides yielding 
the correct bands from a macroscopic approach, our scheme allowed us to classify 
the polarization of each mode as either transverse, longitudinal or mixed. Even 
though the macroscopic approach might fail to yield some modes, namely, those which 
have no coupling to the macroscopic field due to symmetry, they may be recovered 
by extending the notion of macroscopic state, allowing its wavevector to lie beyond 
the first Brillouin zone. The non-locality of the dielectric response may be partially 
accounted for through a magnetic permeability which can then be employed in the 
calculation of the optical properties of the system. We compared the band structure 
obtained through this local approximation to the exact results and we obtained partial 
agreement close to T for those bands that do couple to long-wavelentgh fields. In 
particular, we showed that this approach can yield a negative dispersion in frequency 
regions where both the permitivity and permeability are negative. Nevertheless, we 
discussed some shortcommings of the local approach. Although here we presented 
results for dispersionless transparent dielectrics, our calculation does not require that 
the materials that make up the system be non-dispersive nor non-dissipative, so that 
calculations for real dielectrics and metals may be performed with the same low 
computational costs. 
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